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Abstract 

Linear stability of multi-vector-soliton bound states in the cou- 
pled nonlinear Schrodinger equations is analyzed using a new tail- 
matching method. Under the condition that individual vector soli- 
tons in the bound states are wave-and-daughter-waves and widely sep- 
arated, small eigenvalues of these bound states that bifurcate from the 
zero eigenvalue of single vector solitons are calculated explicitly. It 
is found that unstable eigenvalues from phase-mode bifurcations al- 
ways exist, thus the bound states are always linearly unstable. This 
tail-matching method is simple, but it gives identical results as the 
Karpman-Solev'ev-Gorshkov-Ostrovsky method. 
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1 Introduction 

Nonlinear optics and fiber communication systems are advancing very rapidly 
these days. In this process, a widely-used mathematical model is the coupled 
nonlinear Schrodinger (NLS) equations which govern pulse propagation in 
birefringent fibers [1, 2, 3]. Similar equations with a saturable nonlinearity 
also govern the interaction of two incoherently-coupled laser beams [4, 5]. 

"This work was supported in part by a NASA EPSCoR grant. 



These equations arise in water-wave interactions as well [6, 7]. Solution prop- 
erties of the coupled NLS equations have been examined extensively in the 
past ten years. It is known that these equations admit single- hump vector 
solitons due to a nonlinear mutual trapping effect [8, 9]. When these soli- 
tons are perturbed, they undergo long-lived internal oscillations [10, 11, 12]. 
When they collide with each other, a fractal structure can arise in the pa- 
rameter space [13, 14]. 

Multi-vector-soliton bound states also exist in the coupled NLS equations 
[15]. These states are pieced together by several single- hump vector soli- 
tons. They received much attention because of several reasons. First, in 
fiber communication systems, pulse-pulse interference impairs the system 
performance. If multi-soliton bound states exist, these solutions would have 
implications to system designs. Second, the existence of such bound states 
is noteworthy because they can not exist in the NLS equation [16]. Thirdly, 
these states are closely related to similar states in incoherent laser beams, 
which have been observed experimentally [17]. 

After the numerical discovery of these multi-vector-soliton bound states in 
[15], their analytical construction was made in [18] by a tail-matching tech- 
nique. Their linear-stability problem was studied later in [21] by an extension 
of the Karpman-Solev'ev-Gorshkov-Ostrovsky (KSGO) method [19, 20], and 
small eigenvalues bifurcated from the zero eigenvalue of single vector solitons 
were calculated. These calculations show that multi-soliton bound states are 
always linearly unstable. But this KSGO method is quite involved, thus a 
simpler technique for the calculation of eigenvalue bifurcations is called upon. 

In this paper, we use a new tail-matching method to analyze the linear 
stability of two-vector-soliton bound states in the coupled NLS equations. 
Under the condition that individual vector solitons in these bound states are 
wave-and-daughter- waves (i.e., one component is much larger than the other 
one), and are widely separated, small eigenvalues of these bound states that 
bifurcate from the zero eigenvalue of single solitons are calculated. These 
small eigenvalues are all the non-zero discrete eigenvalues of the two-soliton 
bound states. We found that unstable eigenvalues from phase-mode bifurca- 
tions always exist, thus the bound states are always linearly unstable. The 
present technique is much simpler, but it gives identical results as the KSGO 
method [21]. 
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2 Two-vector-soliton bound states: a review 



The coupled NLS equations 

iA t + A xx + (\A\ 2 + (3\B\ 2 )A = 0, (2.1) 

iB t + B xx + (\B\ 2 + [3\A\ 2 )B = 0, (2.2) 

govern optical pulse propagation in birefringent fibers [1, 2, 3]. Here (3 is the 
cross-phase- modulational coefficient. When = or 1, these equations are 
integrable by the inverse scattering method [16, 22]. 

These equations admit solitary-wave solutions of the following form: 

A(x, t) = r{x)e^ 2t , B{x, t) = i?(x)e^\ (2.3) 

where co and O are frequencies, and the real-valued amplitude functions r(x) 
and R(x) satisfy the ordinary differential equations (ODEs): 

r xx - u?r + (r 2 + l3R 2 )r = 0, (2.4) 

R xx -n 2 R+ (R 2 + (3r 2 )R = 0. (2.5) 

After a simple rescaling of variables, we normalize u = 1. Since Eqs. (2.1) 
and (2.2) are Galilean-invariant, moving solitary waves can be readily de- 
duced from the stationary ones (2.3) (see [23]). 

Solitary waves in Eqs. (2.4) and (2.5) have been studied extensively before 
(see [9] and the references therein) . It has been shown that for any frequency 
SI € [fi c , l/fi c ] where 

n^VTTW-i^ (26) 

this ODE system admits a unique, single-hump, and positive vector-soliton 
solution which is symmetric in both r and R components. We call this 
solution the fundamental vector soliton, and denote it as [ro(x), Rq(x)]. The 
asymptotic behavior of this fundamental soliton at infinity is 

r (x) ^ce-W, R (x) — ► Ce~ n ^, \x\ oo, (2.7) 

where c and C are tail coefficients. When £1 is close to the boundary value Cl c , 
Rq(x) <C ro(x), c 2\/2, and C <C 1; if £1 is close to l/f2 c , ro(x) <C i?o(x), 
c <C 1, and C « 2y/2/Q c . These vector solitons with i?o *C ro or ro <C Rq 
are the so-called wave-and-daughter-waves. 

Two-vector-soliton bound states also exist in the ODE system (2.4) and (2.5) 
[9, 15, 18, 21]. These bound states look like two single-humped vector solitons 
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glued together, while the two solitons are in-phase in one component, and 
out-of-phase in the other component. The in-phase component of the bound 
states are symmetric around the bound-state center, and the out-of-phase 
component are anti-symmetric around the bound-state center. In the limit 
of large soliton separation, these bound states approach a superposition of 
two fundamental solitons (to the leading order): 

r(x) — ► r (x) =F r (x — A), (2.8) 

R(x) — ► Rq(x) ± Ro(x - A), (2.9) 

where the separation A S> 1. These widely-separated bound states exist 
in two parameter-regions [9, 39]: (i) « U c or l/fi c , and < (3 < 1; (ii) 
Q. ~ 1, and f3 > 0. In the first region, the bound states look like two wave- 
and-daughter-waves glued together; while in the second region, the bound 
states look like two nearly-equal-amplitude vector solitons glued together. 
In this article, we only consider the bound states in the first region. In this 
region, the spacing between the two wave-and-daughter-waves in the bound 
state is given by the formula [18, 21] 

. lnc 2 -ln0 2 C 2 

A= i-n ■ (210> 

Below, the bound states with the minus sign in (2.8) and the plus sign in (2.9) 
will be termed type-I bound states, while the ones with the plus sign in (2.8) 
and the minus sign in (2.9) termed type-II bound states (as we have done 
in [21]). These bound states at parameter values (3 = 2/3 and f2 = 0.85 are 
displayed in Fig. 1. A comparison between the analytical spacing formula 
(2.10) and numerical values has been made in [21], and excellent agreement 
has been obtained. 



Analytical construction of multi-soliton bound states in a general nonlinear 
wave system was made in [18] using a tail-matching method, and the spacing 
formula (2.10) for the coupled NLS equations was derived only as a special 
case (see [24] for an application of this method for the construction of other 
types of multi-pulse bound states). Below, we use a simplified version of 
[18] 's method to construct multi-vector-soliton bound states in the coupled 
NLS equations [i.e., (2.4) and (2.5)], and reproduce the spacing formula 
(2.10). There are two reasons for our doing this: (i) to highlight the key 
ideas in the tail-matching method for the construction of multi-pulse bound 
states; (ii) to motivate a similar tail-matching idea for the linear-stability 
analysis of multi-pulse bound states (see Sec. 3). 

Suppose we have a bound state of two vector solitons located at x = and 
A (>• 1). As A — > oo, the leading-order asymptotics of this bound state is 

r(x) — ► r (x) + sir (x - A), R(x) — ► R (x) + s 2 R (x - A), (2.11) 



4 




Figure 1: Numerically obtained stationary two-vector-soliton bound states 
at (3 = 2/3 and Q = 0.85: (a) type-I state; (b) type-II state. Analytical 
approximations by Eqs. (2.8), (2.9), and (2.10) are indistinguishable from 
the numerical curves and thus not shown. 



where si and s 2 are sign-constants and are either 1 or —1. Our task is to 
determine values of s\, s 2 and the spacing A. For this purpose, we consider 
the bound state in two x-regions: — oo < x <C A, and <C x < oo. Since 
the treatments for these two regions are the same, we only look at the first 
region — oo < x <C A. In this region, the bound state is a slightly perturbed 
fundamental vector soliton, i.e., 



r(x) = r (x) + r(x), R(x) = Rq(x) + R(x), 



(2.12) 



where r, R <C 1. The actual forms and sizes of f and R are not important 
in this analysis, but their asymptotics in the region x ~ ^A » 1 is crucial. 
This asymptotics can be obtained by matching [r(x), i?(x)]'s expressions 
(2.12) with their asymptotics (2.11). This is the key idea of the method. 
This matching gives the leading-order asymptotics of (f, R) as 



f{x) 




R(x) 





sxroix - A) 
s 2 R (x - A) 



s\ce 
S2Ce n(x-A) 



x ~ -A » 1. 



As x —>■ — oo, (f, R) — > 0. 



(2.13) 



When Eq. (2.12) is substituted into (2.4) and (2.5), and higher-order terms 
dropped, the linearized equations for perturbations (f, R) are found to be 



L 



r 
R 



= 0, 



(2.14) 
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where the linearization operator L is 
L = 



d xx -l + 3r 2 (x) + PR 2 (x) 2r (x)R (x) 

2r (x)R (x) d xx -n 2 + 3R 2 (x) + (3r 2 (x) 



, (2.15) 



which is self-adjoint. Eq. (2.14) has one localized homogeneous solution 
[r' (x), R' (x)] T when /3 / and 1. Here the prime is the derivative, and the 
superscript "T" is the transpose of a vector. Thus, in order for the solution 
(r,R) of Eq. (2.14) to exist, a solvability condition 



[r' (x),R' (x)]L 

-oo 



r 
R 



dx = 



(2.16) 



must be satisfied. Here x m ~ |A. This condition can be simplified through 
integration by parts as 



r'(x)r' (x) - r{x)r'^{x) + R'(x)R' {x) - R(x)R'^x) 



0. 



(2.17) 



Finally, substitution of the asymptotics (2.7) and (2.13) into the above con- 
dition gives 

(2.18) 



3 -(i-n)A = _; 



This equation readily shows that in order for the bound state to exist, we 
must have s\S2 = — 1. In that case, the spacing formula (2.18) becomes 
exactly the same as (2.10). Hence the above simplified tail-matching method 
reproduces the results from the more general analysis in [18]. 

The relative errors of the leading-order bound states (2.8), (2.9), and the 
spacing formula (2.10) have been discussed in [18], and these errors are 
0(e _A , e~ nA ). For the bound states, we have 



-A -HA> 



r(x) = [r (x) T r (x - A)] [l + 0(e~ A , e 
R(x) = [R (x) ± R {x - A)] [l + 0{e- A , e 



!1A^ 



(2.19) 
(2.20) 



These error estimates can also be obtained as follows. Let us consider the 
type-I bound state. Write 

r(x) = tq{x) — ro(x — A) + f(x), (2-21) 

R(x) = R (x) + R {x - A) + R(x), (2.22) 

where r,R <C 1. Substituting these equations into (2.4) and (2.5), we find 
that in the region — oo < x <C A, 

[3rg(x) + pR 2 (x)] r (x - A) + 2(3r (x)R (x)R (x - A) 
[3R 2 (x) + (3r 2 (x)] R (x - A) - 2(3r (x)R (x)r (x - A), 

(2.23) 



f 




R 
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Here terms which are higher-order than the ones kept in (2.23) have been 
dropped. A similar equation can be obtained in the region <C x < oo. 
Inspection of these equations shows that for all values of x, (f , R) is expo- 
nentially small compared to the leading-order terms in (2.8) and (2.9), and 
the relative errors are 0(e' A ,e- nA ) as shown in (2.19) and (2.20). 

The spacing formula (2.10) together with the results in [9] reveals that, 
in order for A > 1, we must have < (3 < 1. In addition, we must 
have either C/c < 1, fl < 1, or C/c > 1, O > 1. In the former limit, 
Rq(x) <C vq(x), and ps f2 c ; while in the latter limit, vq(x) <C Rq(x), and 
ps l/Q c . In both limits, the bound states look like two wave-and-daughter- 
waves glued together. These two limits are actually equivalent through a 
scaling of variables (the so-called reciprocal relation in [9]). Thus in the 
rest of this article, we will just consider the former limit where 17 ps Cl c . 
Fundamental solitons in this limit have been perturbatively determined in 
[9]. Utilizing those results, we readily find that the asymptotic formula for 
the spacing A is 

A^y-^-{-]n(n-fi c )+tf + o(l)}, n^n c , (2.24) 
where the constant K is 



K = (3- 20 c )ln2- 21nO c + hi7, 



and 
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[1 - n 3 c ) J^sech^xdx 



2Q c J™ OD sech 2nc xdx 



(2.25) 



3 Linear-stability analysis of two-vector-soliton 
bound states 



To study the linear stability of the above two-vector-soliton bound states, 
we perturb these states as 

A = e u {r(x) + Mx)e iXt + ^ 2 *(x) e - iA *'} , (3.1) 

B = e in ' 2t {r(x) + Mx)e iXt + ^t{x)e~ iXH ] , (3.2) 

where [r(x), R(x)] is a two-vector-soliton bound state, ipk (1 < k < 4) are 
infinitesimal disturbances, A is the eigenvalue, and ^* is the complex con- 
jugate of tp. Substituting these equations into (2.1) and (2.2) and dropping 
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higher-order terms, we get the following eigenvalue relation 



= A*, 



(3.3) 



where the linearization operator C is 




1 + 2r 2 + fJR 2 



-prR 
a xI - Q 2 + 2R 2 + f3r 2 
-R 2 




llrR 
-/3rR 



-dxx + 1 - 2r 2 - f3R 2 
/3rR 
-0rR 



■? - fir 2 , 
(3.4) 



and ^ = (ipi, ip2, "03, ^4) T - The spectrum of C contains all information on 
the linear stability of the two-vector-soliton bound state. If C has eigenvalues 
with Im(A) < 0, then the bound state is linearly unstable; otherwise, it is 
linearly stable. Obviously, the continuous spectrum of C always lies on the 
real axis, thus we only need to look at the discrete eigenvalues of C. Notice 
that C 2 is self-adjoint, thus £'s eigenvalues are either purely real, or purely 
imaginary. In addition, if A is an eigenvalue of £, so are —A. Hence Cs 
eigenvalues always come in pairs on the real or imaginary axis. 

It is insightful to view £'s discrete eigenvalues for the bound states (r, R) 
in the perspective of eigenvalue bifurcations from a similar operator Cq for 
fundamental vector solitons (ro,i?o)- Here 



The spectrum structure of Cq has been determined completely in [12, 25]. 
This operator has 6 or 8 discrete eigenvalues (multiplicity counted), depend- 
ing on whether an internal mode exists or not. The zero eigenvalue always 
has multiplicity 6 — three from position and phase invariances (Goldstein 
modes), and the other three from velocity and frequency (or equivalently, 
amplitude) variations. When an internal mode exists, a pair of real eigen- 
values of opposite sign are present as well. If two vector solitons form a 
widely-separated stationary bound state (A S> 1), both solitons must be 
either wave-and-daughter- waves (with Q Q. c and < (3 < 1), or having 
nearly equal amplitudes in the two components (with Q 1 and (3 > 0) 
[9, 39] (only the former type of bound states is studied in this paper). It 
has been shown in [25] that wave-and-daughter-waves do not admit internal 
modes. For instance, single vector solitons in Fig. l's bound states do not 
have internal modes [25] . Thus Cq has only a single discrete eigenvalue zero 
of multiplicity 6. In a bound state of two wave-and-daughter-waves, C will 
then have 12 discrete eigenvalues (multiplicities counted) — double that for 
a single wave-and-daughter-wave. Here, no new discrete eigenvalues of C 
can be generated from edge bifurcations of Cq because the edge points of 
£o' s continuous spectrum are not in the continuous spectrum. Now the zero 
eigenvalue of C still has multiplicity 6. Another six eigenvalues of C must bi- 
furcate from zero due to tail interactions between solitons. Tracking of these 



£o - C\{ r =r ,R=R )- 



(3.5) 
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six bifurcated eigenvalues will then provide a complete characterization of 
linear stability for these two-soliton bound states. 

In this section, we propose a new and general tail-matching method to derive 
explicit analytical formulas for the six eigenvalues of C that bifurcate from 
zero, under the condition that the individual vector solitons in the bound 
state are widely separated. The idea of this method is to perturbatively de- 
termine the bifurcated eigenstate around each vector soliton. By matching 
their tail asymptotics in the center region of the eigenstate, their asymp- 
totics in that region can then be explicitly obtained. Finally by utilizing the 
solvability conditions for the eigenstate, analytical formulas for the bifur- 
cated eigenvalues will be derived. These formulas turn out to be exactly the 
same as those obtained by the KSGO method [21], but the present method 
is much simpler. 

The bifurcated eigenvalues are of two types. One type consists of a pair 
of eigenvalues which bifurcate from the position-related zero eigenvalue. At 
infinite soliton separation, these eigenstates are equal to the sum of two 
position-induced Goldstein modes separated infinitely apart. The other type 
consists of two pairs of eigenvalues which bifurcate from the phase-related 
zero eigenvalue. At infinite soliton separation, these eigenstates are equal to 
the sum of two phase-induced Goldstein modes separated infinitely apart. 
These two types of eigenvalues have their counterparts in the linear stability 
analysis by the KSGO method [21]. 

Below, we consider eigenvalue bifurcations in type-I and II bound states. It 
turns out that eigenvalues of type-II bound states are simply equal to those 
of type-I states multiplied by i (see also [21]). Thus calculations for only 
type-I states will be presented. These states are anti-symmetric in r, and 
symmetric in R, around the center of the bound states [see Fig. 1(a)]. When 
A -»■ oo, 

r(x) — ► r (x) — r (x — A), R(x) — ► Ro(x) + R (x — A). (3.6) 

Bifurcations of position- and phase-related eigenvalues are studied separately 
next. 

3.1 Bifurcation of position-related eigenvalues 

When A — > oo, the eigenstate bifurcated from the position-related zero 
eigenvalue is a sum of two position-induced Goldstein modes of fundamental 
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vector solitons: 



¥(ar) — >^o(x) + ^o(x), A ->■ oo, 



(3.7) 



where 



*o(x) 



r&(x) 
r&(x) 
i^(z) 



r' (x-A) 
r' Q (x — A) 
-R> (x-A) 
-R' (x - A) 



(3.8) 



The first two components of ^ are anti-symmetric around x = |A, and the 
last two components are symmetric around x = ^A. Note that in the same 
limit, the eigenstate ^o(x) — ^o(x) is simply [r'(x), r'(x), R'(x), R'(x)] T , 
which is the un-bifurcated Goldstein eigenmode with eigenvalue zero and is 
thus not considered. 



In the limit A » 1, we consider the bifurcated eigenstate in the region 
— oo < x <C A, and expand it as a perturbation series in powers of the small 
eigenvalue A: 

^(x) = y (x) + \y 1 (x) + \ 2 y 2 {x) + o(\ 2 ), -oo<x«A. (3.9) 

In this region, we also expand 

£ = £ + A 2 £ 2 + o(A 2 ). (3.10) 

It is noted that when A>1, single solitons in the bound state considered are 
wave-and-daughter-waves whose two components have different orders (one 
component asymptotically much smaller than the other). This fact certainly 
has implications in the stability analysis. In particular, different components 
of <I> and C in the region — oo < x <C A may have slightly different orders of 
magnitude. Thus, a perturbation expansion with a more-detailed ordering 
of different components than that in (3.9) and (3.10) might be needed. But 
as we will see next, (3.9) and (3.10) still work. 

Now we substitute the perturbation expansions (3.9) and (3.10) into the 
eigenvalue equation (3.3). At O(l), we get Cq^q{x) = which is satisfied 
automatically. At O(A), we get 

£ *i = *o, (3.H) 

whose solution is 

^i(x) = [^xr (x), ~xr {x), ^xRq(x), -^xR (x)] T . (3.12) 

Note that this function is an inhomogeneous solution of Eq. (3.11). In 
general, ^\ should also include the homogeneous solutions which are a linear 
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combination of the three Goldstein modes: *&o(x) in (3.8), and [$q^(x), 

(2) 

$q (x)\ in (3.26). But the ^o(x) term in can be combined with the O(l) 
term in the expansion (3.9), and the other two Goldstein modes 

$q 2 ^(x)] are phase-related and do not arise here. Thus the solution of Eq. 
(3.11) can be taken as (3.12) without loss of generality. 

At <3(A 2 ), we get 

£0^2 = *i - £ 2 1 , o- (3.13) 

It is more convenient to express £2^0 in a different form. Recall that the 
un-bifurcated position-related Goldstein mode ^f g (x) of £ has the leading- 
order asymptotics ^o(x) — ^q(x) for A> 1. When this asymptotics and £'s 
expansion (3.10) are substituted into the Goldstein-mode relation g (x) = 
0, we find that asymptotically, 

A 2 £ 2 ^o = £o^o- (3.14) 



Thus Eq. (3.13) can be rewritten as 



£ 



^2 + A" 2 ^ = *i- (3.15) 



Note that £0^0 is 0(X 2 ), thus the above equation is not dis-ordered. 



The solvability condition of Eq. (3.15) will produce formulas for the eigen- 
value A. To do this, we need the asymptotics of function ^2 in the region 
x ~ > 1, which we derive using the tail- matching idea. We have known 
that 'f's leading-order asymptotics at A ^ 1 is given by (3.7) for all values 
of x. Combining this asymptotics with the perturbation expansion (3.9) and 
solutions (3.8) and (3.12), we see that 

*2(*) ^o(x), x ~ ^A > 1. (3.16) 
Of course, ^2 0*0 — ► when x — > —00. 

The homogeneous equation of (3.15) has three linearly independent solutions 
which are the Goldstein modes ^o(x) in (3.8) and [$o^(x), 3>q 2 ^(x)] in (3.26). 
Thus Eq. (3.15) has three solvability conditions. These conditions can be 
readily derived by noting that diag(l, — 1, 1, — l)£o is self-adjoint. It turns 
out that two of the solvability conditions induced by the Goldstein modes 
[^o^j^o^] are satisfied automatically. The remaining solvability condition 
reads, 




(*o|diag(l, -1, 1, -1)|£ (*2 + \- 2 * ))dx 



<*o|diag(l,-l,l,-l)|*i)dx, (3.17) 

-00 
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where x m ~ ^A, (-| and |-) are the Dirac ket and bra notations [26]. Inte- 
grating by parts to its left-hand-side and simplifying its right-hand-side, we 

get 

(^oK^ + A-^o*))- (^ox|(^2 + A- 2 l-o)) m =--/ (r 2 + Rl)dx. 

(3.18) 

The left-hand-side of (3.18) can be calculated using the asymptotics (2.7), 
(3.16) and equation (3.8), while the integral on the right-hand-side of (3.18) 
is asymptotically equal to a similar integral but with the upper limit x m 
replaced by oo. After these calculations, the eigenvalue A is finally found to 
be 

M + N x ' 

where M and N are the masses of the r and R components in a fundamental 
vector soliton: 

/OO POO 
rgdx, N= Rldx. (3.20) 
-oo J —oo 

We immediately see that formula (3.19) is identical to the one derived in [21] 
using the KSGO method. Thus, the present tail-matching method has the 
same accuracy as the KSGO method, but is only simpler. Recall that we 
only consider the « Q c (< 1) limit. In this limit, the type-I state flips sign 
in its larger component, and does not flip sign in its smaller component [see 
Fig. 1(a)]. According to formula (3.19), this position-related eigenvalue A is 
real, thus it does not create instability. A comparison between the analytical 
formula (3.19) and numerical values at (5 = | has been made in [21], and 
excellent agreement has been obtained. 

Formula (3.19) shows that A = 0(e~5 A ). When fi-»fi c andO</?<l, 
A — > oo (see Sec. 2). In this limit, the asymptotic formula for A can be 
obtained more explicitly from (2.24) and (3.19) as 

a 2 — >a(n-n c )^, n^n c , (3.21) 

where the constant a is 

Here Q c and 7 are defined in Eqs. (2.6) and (2.25). 
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3.2 Bifurcation of phase-related eigenvalues 



Calculations for the bifurcation of phase-related eigenvalues are quite similar 
to that done above. In this case, the eigenmode has the following asymp- 
totics: 

9(x) — >& (x) + &o(x), A 



oo, 



where 



$o{x) = ^q\x)+S^q\ 



X 



Mx) = n'(x-A)-6^'(x-A), 

r (x) 
-r (x) 








*£ 2) (*) 






Ro(x) 
-Rq(x) 



(3.23) 

(3.24) 
(3.25) 

(3.26) 



and 5 is some constant (to be determined). The first two components of 
this mode are symmetric around x = |A, and the last two components are 
anti-symmetric around x = |A. Note that the eigenstate with asymptotics 
$>o(x) — ^o(^) is simply [r(x), —r(x), SR(x), — 5R(x)] T , which is the un- 
bifurcated Goldstein eigenmode with eigenvalue zero and is not a concern. 



Next, we construct a perturbation-series solution for \l/(x) in the region 
— oo < x <C A. The perturbation series is similar to (3.9): 



= $ (x) + A*i(x) + \ 2 <S> 2 (x) + o(A 2 ), 



-oo < x < A, (3.27) 



where &o(x) is given by (3.24). Substituting this expansion and (3.10) into 
the eigenvalue relation (3.3), the 0(1) equation is satisfied automatically. At 
0(A), we get 

A)$i = $o, (3-28) 



whose solution is 



$i(a:) 



ld_ 



r (x;uj,n) 
r (x;uj,n) 
Ro(x; uj, Q) 
Ro(x; (jj, f2) 



+ 



_5_d_ 



U>=1 



ro(x;w,$7) 
r (x;tv,£l) 
Rq(x; to, fi) 
i?o(a^; $7) 



(3.29) 



Here (ro,-Ro) is the fundamental vector soliton obtained from ODEs (2.4) 
and (2.5) without rescaling of u> = 1. Again, the homogeneous solution of 
Eq. (3.28), which is a linear combination of Goldstein modes ^o(x) in (3.8) 



and [$>£\x),&q>(x)} in (3.26), is not included because the latter modes can 
be absorbed into the 0(1) term in the perturbation expansion (3.27), and 
the former mode does not arise. 



,(2) 
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At <3(A 2 ), we get 

£ $ 2 = - £ 2 $ . (3.30) 

Again, utilizing the un-bifurcated phase-related Goldstein modes of L with 
asymptotics <&o(x) — &o(%) — similar to what we have done in Sec. 3.1, we 
can rewrite £2^0 so that Eq. (3.30) becomes 



$2 + A 2 l>o 



$1. (3.31) 



This equation has three solvability conditions induced by the three Goldstein 
modes in the homogeneous solution. The condition induced by the mode 
^o(x) in (3.8) is satisfied automatically. The other two conditions are 

f'V^ldiagU, -1, 1, -l)|£o($2 + \- 2 $o))dx 
J — 00 

= f m <*( fc) |diag(l,-l, l,-l)|$i)dx, (3.32) 

J — 00 

where x m ~ ^A, and k = 1,2. Integration by parts simplifies these condi- 
tions as 

"(^K^ + A- 2 ^,)) - <^S?|(*2 + A- 2 $o))]^ = + ^M Q , 

(3.33) 

and 

1 S 

'(^o ] \^2 X + A- 2 !^)) - ($S|($2 + A- 2 <&o)>]^ = + (3.34) 

To calculate the left-hand-sides of the above two equations, we need the 
asymptotics of function $2(2^) m the region x ~ ^A. This asymptotics can 
be obtained by comparing ^(x)'s asymptotics (3.23) with its perturbation 
expansion (3.27). This comparison shows that ^(x) must have the asymp- 
totics 

$ 2 (x) - ^ (x), x ~ ^A > 1. (3.35) 

When this asymptotics as well as (2.7) is substituted into Eqs. (3.33) and 
(3.34) and parameter 5 eliminated, the eigenvalue A is then given by the 
quartic equation 

4 Wc\N n - M^e- A 2 W^e^ 

M^N n - MqN^ M u N n - MqN^ ' 1 ' ' 

Again, this formula for phase-related eigenvalues is identical to that ob- 
tained in [21] by the KSGO method. As pointed in [21], this formula shows 
that a two-vector-soliton bound state always has one phase-related unsta- 
ble eigenvalue which bifurcates from zero, thus is always linearly unstable. 
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Comparison between this formula and numerical values has also been made 
in [21], and excellent agreement has been observed. 

Formula (3.36) shows that phase-related eigenvalues A = 0(e~2 A ) 5 the same 
as position-related eigenvalues [see Eq. (3.19)]. The asymptotics of these 
eigenvalues in the limit f2 — > il c can also be obtained from (2.24) and (3.36), 
but is not pursued here. 

Next, we briefly discuss the linear stability of type-II vector-soliton bound 
states [see Fig. 1(b)]. In the limit of large separation, these solitons have 
the asymptotics 

r{x) — ► r (x) + r (x - A), R(x) — > R (x) - R (x - A) . (3.37) 

Repeating the above analytical calculations, we find that the eigenvalues 
for type-II states are equal to those of type-I states multiplied by i. Thus 
type-II states are always linearly unstable as well. But different from type-I 
states, the instability of type-II states in the limit O « fi c (< 1) is caused 
by two unstable eigenvalues: one related to position-mode bifurcations [see 
(3.19)], and the other one related to phase-mode bifurcations [see (3.36)]. 
This result agrees with that by the KSGO method as well as the numerical 
result [21]. 



4 Discussion 



In this paper, we have analytically studied the linear stability of two- vector- 
soliton bound states in the coupled NLS equations by a new tail-matching 
method. Under the condition that the two vector solitons are wave-and- 
daughter-waves and are widely separated, we have calculated small eigenval- 
ues of these bound states that bifurcate from the zero eigenvalue of single 
vector solitons. These small eigenvalues calculated are all the discrete non- 
zero eigenvalues of the bound states. We have shown that these bound states 
are always linearly unstable due to the existence of one unstable phase- 
induced eigenvalue. The analytical formulas for eigenvalues derived from 
this tail-matching method turn out to be exactly the same as those from the 
KSGO method [21], but the present method is much simpler. Even though 
our calculations were performed for two-vector-soliton bound states, these 
calculations can apparently be generalized to n-vector-soliton bound states. 

This tail-matching method for the linear stability analysis of multi-soliton 
bound states is apparently quite general, and it can be applied to a wide 
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range of other wave systems where similar multi-pulse bound states have 
been reported [18, 27, 28]. In addition, this method should be applicable 
to the stability analysis of other types of multi-pulses from nonlocal bi- 
furcations in coupled-NLS-type equations [24]. Preliminary linear stability 
results through numerical studies shows that such multi-pulses are also lin- 
early unstable [24], consistent with the results of the present article. This 
tail-matching method can also be used to calculate eigenvalue bifurcations of 
multi-pulses from internal modes (non-zero eigenvalues) of single pulses — a 
bifurcation which incidentally does not arise in the present problem because 
internal modes do not exist in single wave-and-daughter- waves [25]. 

Under what conditions can this tail-matching method give useful results? 
The main condition is that the individual pulses in the multi-pulse bound 
state are widely separated. This condition will dictate in what parameter 
regions such multi-pulses can exist, and thus tail-matching can proceed. For 
instance, for the coupled NLS equations (2.1) and (2.2), the spacing formulas 
(2.10) and (2.24) dictate that widely-separated multi-vector-soliton bound 
states (A » 1) exist when O is near the local bifurcation boundary f2 = Cl c . 
This is precisely the parameter region where our tail-matching linear stability 
analysis is performed. 

The two-vector-soliton bound states studied in this paper reside outside the 
continuous spectrum of the corresponding linear-wave system. It is known 
that multi-pulse embedded solitons residing inside the continuous spectrum 
exist in various wave systems as well [29, 30, 31, 32, 33]. An interesting open 
issue is whether the tail-matching method can also be applied to the linear 
stability analysis of such multi-pulse embedded solitons. In the third-order 
NLS equation, it has been shown numerically that such embedded solitons 
are all linearly stable, but nonlinearly semi-stable [33]. 

Lastly, we relate this tail-matching method to other existing techniques for 
the linear stability analysis of multi-pulse bound states. Currently, the 
following techniques exist: the KSGO method [21], the dynamical-system 
method [34, 35], and the effective-interaction-potential method [36, 37, 38]. 
The present method is asymptotically accurate. It gives the same results 
as the KSGO method, but is much simpler. The dynamical-system method 
can count the number of unstable eigenvalues, or express the eigenvalues 
as the zeros of the Evans function. But it generally does not produce ex- 
plicit formulas for eigenvalues. The effective-potential method can only cap- 
ture position-related eigenvalue bifurcations, not phase-related eigenvalue 
bifurcations. (For the coupled NLS equations, phase-related eigenvalues are 
more important). In view of this comparison, we feel that the tail-matching 
method for the linear stability of multi-soliton bound states is very promis- 
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mg. 
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